Auxiliary— field Monte Carlo for Quantum Spin and Boson Systems 
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We describe a new algorithm for the numerical simulation of quantum spin and boson systems. The 
method is based on the Trotter decomposition in imaginary time and a decoupling by auxiliary 
Ising spins. It can be applied, in principle, to arbitrary (random) spin systems, however in general 
it suffers from the "minus-sign problem" . This problem is absent in the case of the Ising model in a 
transverse field in arbitrary dimensions and geometries. We show test results for the spin-1/2 XY 
model, the one-dimensional transverse Ising model with disorder, and the phase transition induced 
by a transverse field in the two-dimensional Ising model. 
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I. INTRODUCTION 

Quantum Monte Carlo (QMC) methods have con- 
tributed to much of our recent knowledge of the prop- 
erties of interacting quantum mechanical spin systems, 
and closely related hard-core boson models. Green Func- 
tion Monte Carlo (GFMC) is one powerful class of ap- 
proaches which project out the lowest energy many-body 
eigenfunction. World-line (WLQMC) algorithms consti- 
tute another generic class of QMC approaches, and allow 
the evaluation of finite temperature properties. Recently, 
very substantial improvements to WLQMC, the loop ||] 
and continuous time techniques, have been developed. 

Both GFMC and WLQMC most commonly use a basis 
labeled by the boson occupation number or position, or, 
analogously, the z component of spin, in space and imag- 
inary time. The key feature of the approaches is that 
eigenvalues of the original operators in the Hamiltonian 
describe the Monte Carlo configurations. In contrast, the 
preferred techniques for QMC simulations of interacting 
fermions || , involve the introduction of an auxiliary field. 
The original fermion operators are integrated out, and 
the simulation takes place in the space of this abstract 
auxiliary field. 

In this paper we will introduce a new auxiliary field 
QMC method for interacting quantum mechanical spins 
and boson systems. Why is such an algorithm interest- 
ing? WLQMC and GFMC approaches have very sig- 
nificant weaknesses, including extremely long correlation 
times and restrictions on the observables which can be 
measured. While loop algorithms Q have addressed this 
issue, their efficiency remains problematic in several im- 
portant cases, for example when interactions are longer 
range, or disorder is present. Therefore, continued algo- 
rithm development is desirable. 

The organization of this paper is as follows: We will 
first introduce the Ising model in a transverse field, and 
briefly review the key issues in its properties. We will 
then describe how an auxiliary field algorithm can be 
constructed for this model. Although related conceptu- 
ally to fermion QMC, it differs considerably from analo- 



gous fermion techniques in that the resulting traces are 
over independent single site problems, avoiding the neces- 
sity to evaluate the determinants of large matrices in the 
fermion case. We then give results of our simulations, in- 
cluding a comparison of the approach with existing tech- 
niques. We conclude by describing another interacting 
spin/boson model, the boson-Hubbard model. However, 
we show that the sign problem is a serious limitation to 
auxiliary field approaches in this case. 



II. TRANSVERSE FIELD ISING MODEL 

Two quantum spin/hard-core boson problems of con- 
siderable recent issue are the Ising model in a transverse 
field [0, and the boson-Hubbard model ||. The former 
allows one to study in a simple setting many of the key 
qualitative issues in quantum phase transitions in disor- 
dered systems, including the nature of the distribution of 
correlation functions and the shifts in the values of crit- 
ical exponents from the clean limit. The latter offers a 
description of the superconductor-insulator phase tran- 
sition when preformed pairs exist above the transition, 
and in the hard-core limit is also formally identical to the 
quantum mechanical spin-1/2 XXZ Hamiltonian. In this 
section we will describe the Ising model in a transverse 
field, which appears to be the more promising application 
of the auxiliary field approach. 



A. Hamiltonian and Algorithm 



The Transverse Ising model [|],[3| is given by 



h = Justus* - ^ B i s i 



(i) 



Here Sf, a £ {x,y,z}, are the Pauli matrices obeying 
the commutation relations: [S",Sj] — Sij(e a /3 7 S] + S a p). 
The sum runs over pairs of nearest neighbors, (ij). The 
partition function is given by Z — Tr e~^ H . 
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We employ the usual Trotter break-up of non- 
commuting operators in the Hamiltonian: 



(2) 



/V 



n 

i=l \j£NN(i) 
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0(At 2 ) 



with AtL — p. The inner product runs over the Z near- 
est neighbors of i. 

In order to decouple the interaction terms we recall 
that any product of two commuting operators can be 
written as a sum over squares: 



AB = ~ [(A + Bf - (A - B) 2 ] , 



(3) 



and a squared operator can be decoupled by the introduc- 
tion of a Gaussian integration over a classical auxiliary 
field: 



e (A+BY/i = w -i/2 



(4) 



In the present case, one squared operator is sufficient 

SIS* = a [2(P%) 2 - 1] , with P% = {St + aS?)/2 (5) 
and a = ±1. 

Employing that [P^) 2 is a projection operator, 
i.e. {P%) 2k = (P l a J ) 2 , for k = 1,2,---, one immediately 
confirms by Taylor expansion that: 

e- 2QArJ -( p 5) 2 =cosh(2A y ^) 
2 

<7=±1 



(6) 



with cosh(2Ay) = exp(— 2aAr Jy). In order to get real 
variables one chooses a = +1 (—1) for < (> 0). 
Thus a two-valued rather than Gaussian decoupling of 
the interaction is possible, introducing Ising-type auxil- 
iary spins cry = ±1: 



I e -ArJ„ 



E 



(7) 



=±i 



for Jij > 0. This "discrete Hubbard-Stratonovich trans- 
formation" was first introduced by Hirsch in the fermion 
case 0. 

The same decoupling holds for any component of the 
Pauli matrix, hence the extension to XY or (anisotropic) 
Heisenberg models or, equivalently, hard-core bosonic 
models, as will be discussed. 

The decoupling has to be done for each of the L factors 
in (^|) giving the auxiliary spins 0^ (7) an additional index 
I = 1, - ■ ■ , L. As a result, we have a system of non- 
interacting Ising spins in a transverse field Bi coupled to 
an auxiliary longitudinal two-valued field, described by 



(Jij (I), which fluctuates in space and "imaginary time". 
For a given configuration, the original spins are 

trivially described by a direct product of 2 x 2 matrices. 

Note that <Jij(l) is a bond, not site variable, i and j 
denote pairs of sites connected by a Jij ^ 0. Hence, there 
are NZL/2 auxiliary spins where Z is the coordination 
number of the lattice. Even in the classical limit (Bi = 0, 
L=l) the auxiliary spins are not dual to the original ones. 
An exception is the one-dimensional classical case where 
original and auxiliary spins are equivalent. 

As mentioned above, for a given configuration of auxil- 
iary spins the original spins are independent, the Hilbcrt 
space factorizes, and the partition function can be writ- 
ten as: 



N 



Z = E Il Tr *n II e A « ff « Ws . ? e ATBii T| (8) 



{°-y(0} <=1 1=1 V'£NN(») 



E 



2x2 matrix product 
positive definite weight function 



(9) 



That is, Z is now a sum over the NZL/2 auxiliary Ising 
spins with a weight function proportional to the product 
of traces of 2 x 2 matrices, one for each lattice site. 

In the case of the Ising model one can choose all lo- 
cal transverse magnetic field values Bi to be positive or 
change their sign by a local spin rotation, respectively. 
Thus there are only positive matrix elements involved 
and w({(Jij(l)}) can serve as a positive definite weight 
function. This is, however, different, e.g. in the XY model 
where severe sign problems occur even for small systems 
(see Sec. III). 

To simplify the notation for the following, we denote 
every factor in the matrix product (||) by 

Bi(l)= ( e *W)sf j ,±td.s; 

yeNN(i) 

and define the product over I and its cyclic permutations 
as 

Ml) = Bi(l) Bi{l + !)■■■ Bi(L) B t (l) ■ .-BiQ - 1) (11) 

In fact, the values of the traces do not change under 
cyclic permutation in the matrix product, i.e. they do 
not depend on the index I in Ai(l), and we can rewrite 
Z as: 



z = E U^iMi) 



(12) 



The resulting Monte Carlo algorithm is similar to the 
auxiliary-field methods for lattice fermions |j, with the 
key difference that one has the product of traces of 
NZL/2 matrices of dimension 2 to evaluate, instead of 
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the determinant of a single matrix of dimension the spa- 
tial lattice size N. Thus the algorithm scales linearly in 
N. It goes as follows: 

1. One starts at "time slice" 1 = 1, initializes the 
auxiliary field (Tij(l), and calculates the 2x2 ma- 
trix products Ai{l) for each lattice site. The possi- 
ble matrix elements (exp(±Aj,5f ), exp(Ar_B i S'f )) 
needed to determine Bj(Z) in ( |1C| ) are calculated 
once for all at the beginning. 

2. At a fixed time slice I try NZ/2 single-spin flips, 

— > cr^(Z) = —<7ij(l). Such a flip involves the 
two lattice sites i and j, and is accepted according 
to the probability ratio 



P 



Tr Ml) TrA'^l) 
Tr MJ) TrAj(iy 



with the new matrices 

A' i {l)=B'S)B i {l)- x Ml). 



(13) 



(14) 



B[(l) denotes the matrix ( |10| ) with <Jij(l) replaced 
by a^Al). Note that it is not necessary to perform 
the whole product of L matrices but only a few 2x2 
matrix operations for each spin flip. If accepted the 
new matrices A'^l) replace the old ones. 

3. Move to the next time slice, I + 1, by calculating 



A i (l + l)=B i (l)- 1 A i (l)B i (l), 



(15) 



for each lattice site i. 
4. Move to (2). 

After L cycles (2-4) one sweep through the d + 1 dimen- 
sional system is complete. It takes oc NL(Z/2) 2 multi- 
plications. Step 3 leads to round-off errors, in particular 
at large /3, so one has to recompute the matrices Ai(l) 
from scratch from to time, typically after ten time slices. 

We note that the systematic error due to the Trotter 
decomposition can be strongly reduced by a third order 
decoupling. 



3 At(A+B) 



At A/2 AtB AtA/2 

e ' e e ' 



C(Ar 3 



(16) 



While the leading correction in expectation values of her- 
mitian operators is C(Ar 2 ) for both second and third or- 
der decoupling, the prefactors are typically a lot smaller 
in the latter case. The implementation is simple. For- 
mally the matrices (|l(]) are changed to 



BS) ee e A ^r/ 2 | Y[ 



In the product ([11]), however, there are always two of 
such factors, exp(Ari? i S'f /2), adding to exp(ArBiSi), 



and the remaining factor at the beginning of the product 
can be shifted to the end since the trace does not change 
under cyclic permutation. Hence, the Monte Carlo pro- 
cedure remains completely unchanged, and it is sufficient 
to replace each local operator Oi by 

-ArBiSf/2 q gArSiSf/2 



O, 



(17) 



The computational effort for these 2N additional 2x2 
matrix multiplications is negligible. 



B. Observables 

The algorithm allows for the measurement of a variety 
of static and time-dependent observables. Since measure- 
ments for successive time slices are in general correlated 
they are performed after every full sweep over space and 
time. 

Interestingly, in the weak coupling limit (Jy = 0), 
all expectation values become exact, independent on the 
auxiliary field configuration, i.e. without any sampling. 
This is not the case if one samples over the original 
Ising spins. By analogy the auxiliary field approach for 
fermions || exactly solves the noninteracting problem 
without sampling, while world-line approaches |8| do not, 
and still require a full Monte Carlo simulation to get ob- 
servables. 



1. Static correlation functions 

Most static observables can be expressed in terms of lo- 
cal magnetizations and static correlation functions. The 
components of the local magnetization Sf are given by 

< 5 "> = ^ E (il 1 ^ 1 )] Tr {^a)} 

1 V- n ^^^^ Tr{Sf^(l)} 
= — > wilo-iAln) , , \ (18) 

Z s „u TrA(l) V ' 

(.()} 

with the weight function w({aij(l)}) from Eq. (||). That 
means we have to sum up the ratio of traces on the 
r.h.s. of ( p8[ ) over the auxiliary field configurations. Sim- 
ilar equations hold for static correlation functions. In 
short: 



Tr{S?Ml)} \ 
TrA(l) /, 




= ' Tr{SfA(l)}Tr{SfA(l)} ' 
' ' TrA(l)TrA(l) 



TrA(l) 



(19) 

(i^k) (20) 
(21) 



where (■ • -) w stands for the sum over {o~ij(l)} configura- 
tions with proper weight. 
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2. Static Susceptibilities 

Susceptibilities, in general, require the calculation of 
correlation functions in imaginary time. The homoge- 
neous susceptibility for spin component a is defined as 



X a =N 



dr (M a (T)M a (0)) - f3(M c 



with the magnetization operator 



1 N 

M a (r) = e- TH M a e TH , If = E ^ 



(22) 



(23) 



C m {l) = B m {l), C m {l + l) = C m (l)B m {l + l) (29) 

n rn (L) = i, n m (i-i) = B m (i)n m (i). (30) 

For each value of I, the product £ m (l)lZ m (l) = A m (l). 
To check for round-off errors this equality is tested from 
time to time. No significant deviations were found for 
the parameters used. 

Xmn a l so determines the dynamical susceptibility in 
imaginary time from which, in principle, the real time 
dynamics can be extracted by an analytic continuation. 



In our discrete time approach the integral is replaced by 
a sum over time slices, yielding: 

L 

X a = ArN^2(M a (ATl) M a (0)) - (3N(M a ) 2 
i=i 

= ^ E E < 5 ™ W S n(°)) ~ PN(M*)\ (24) 

1=1 n,m 

where S*(l) = S^(AtI) means: 

/l \ /i X" 1 



n 5 ™ n 



(25) 



For the time-dependent correlation function we again 
employ the fact that for a given auxiliary-field configura- 
tion operators for different lattice sites commute and we 
obtain for m =/= n: 



XrnJl) = (5° (Z)5°(0)> 



\ E II TrA(l) Tr[^(/)] Tr [5°(0)] 
Tr [5« (/)] Tr [5»(0)] 



TrA„(l) Tr Ai(l) 



(26) 



where the operators (2x2 matrices) in brackets are de- 
fined as: 

l L 

= l[B m (l')S% l J] B m (0 (27) 

Z' = l l'=J+l 

For the on-site correlation function we obtain similarly: 
'Tr{[S£(Z)] 



Tr„4 m (l) 



(28) 



Eq. (p7p can be written as a product C m {l)S^R. m {l) 
where the matrix products on the left and right hand side 
of the spin operator are calculated iteratively: 
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FIG. 1. Random transverse Ising chain with one disorder 
configuration. 2nd and 3rd order Trotter decomposition. Ex- 
act results (x) from Jordan- Wigncr transformation. 
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C. Results 

1. Random field, random bond transverse Ising chain 

The one-dimensional Ising model in a transverse field 
can be solved exactly using the Jordan- Wigner transfor- 
mation H[hJ. With random bonds and/or random mag- 
netic fields, explicit formulas for finite open chains at 
finite temperatures were given by Young [ Ll[ which can 
be used to test the algorithm. We did simulations of 
open chains of 20 and 100 sites with one disorder con- 
figuration, and calculated energy and nearest neighbor 
correlation functions (see Fig. [l]). The energy values are 
compared with the numerically exact ones. The conver- 
gence with At 2 to the exact values is quite good. The 
third order Trotter break-up leads to significantly smaller 
systematic errors in the energy. In the spin-spin correla- 
tion function, however, the prefactor is somewhat larger. 



2. 2D Transverse Ising Model 

A second application of the algorithm is the phase 
transition induced by a transverse magnetic field in the 
pure (non-random) two-dimensional ferromagnetic Ising 
model. Fig. || shows results for one fixed system size at 
an inverse temperature (3 = 10. All quantities are ex- 
trapolated to At — > 0. 




2.0 3.0 4.0 5.0 



FIG. 2. Longitudinal and transverse magnetization, 
M z and M x , staggered susceptibility, xaf, and energy E 
vs. transverse magnetic field B. 10 x 10 lattice sites, /3 = 10 
with At — > extrapolation. Dashed line: fit M z oc (B c — B) 13 
with B C (T = 0) = 3.08. 

The phase transition is visible in the longitudinal and 
transverse magnetization, M z and M x , as well as in 
the staggered susceptibility, xaf, which shows a kink 
at the transition. The homogeneous susceptibility, how- 
ever, was too strongly fluctuating to give reliable results. 



The energy, E, behaves smoothly at the transition as ex- 
pected for a second order transition. M z vanishes around 
B 3.1. If we assume that the finite system is essen- 
tially in its ground state and take the value of the critical 
magnetic field from high temperature expansions Jl2| , ^3| , 
B c = 3.08, then we can fit M z near the transition by a 
power law and get the exponent (5 = 0.32±0.01. This is in 
remarkably good agreement with the value for the classi- 
cal d+ 1 = 3 dimensional Ising model, (3 = 0.325±0.0015 
|]l4| which should apply at the quantum critical point 
(T = 0). 

Auto-correlation times are short even close to the tran- 
sition point. At B = 3.1 we observe auto-corrclations 
in M z , E, and xaf below 0.1 after one sweep, and in 
M x and \f after four sweeps. The data for Fig. || were 
obtained within about one day of computer time on a 
workstation. 



III. THE BOSON-HUBBARD (XXZ) MODEL 
A. Hamiltonian and Algorithm 



The boson Hubbard model is 

H= -*X>l a i +a l a i) 



+V niUj+U rii(jii — 1) — fi ; 



(31) 



Here a i and a\ are the destruction and creation operators 
for bosons on site i, and n; = a\a i is the number opera- 
tor. The first term is the kinetic energy, and U and V are 
on-site and near-neighbor repulsions between bosons, 
is a chemical potential which controls the density of elec- 
trons per site, (n), on the lattice. 

To briefly illustrate some of the properties of the 
model, consider the ground state phase diagram at V = 
0. As the chemical potential is raised, the density of 
bosons on the lattice increases smoothly from zero. How- 
ever, if U is large, then when the density goes through 
(n) = 1, the chemical potential takes a sudden jump, 
since at that point lattice sites become doubly occupied 
at the cost of the big on-site repulsion U. Similar jumps 
occur at all integer fillings. These jumps represent a gap 
in the many-body energy spectrum, and reflect the fact 
that the system is insulating at strong coupling. If U is 
sufficiently small, the kinetic energy dominates, and the 
gap vanishes. Therefore, as U/t is changed, the boson- 
Hubbard Hamiltonian undergoes a quantum phase tran- 
sition between superfluid and insulating states. Away 
from integer filling, the system is a superfluid at any ra- 
tio of U/t. Nonzero V, or the introduction of disorder, 
likewise allow for interesting new phases at T = 0. 

In the hard-core limit, the occupations are rii = 0, 1. 
Then U drops out of the problem and with the usual 
mappings 
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i\^S+ = (Sf+iS»)/2 
Sr = (S* + iSy)/2 



a i 
m 



Sf + 1/2, (32) 
Eq. (29) transforms into the spin-1/2 XXZ model, 

H = -t/2j2(SfSJ + S^S]) 

m 



-v^{si 

(ij) 



i 



(33) 



An occupied site is hence identified with a spin up, etc., 
and the issues discussed above for the boson-Hubbard 
model can be reformulated in spin language. For exam- 
ple, the competition between superfluid and insulating 
phases corresponds to the between magnetic order in the 
XY and Z directions. 

As mentioned above, the auxiliary-field decoupling (0) 
can be performed for different spin components indepen- 
dently. In the case of the y component it is useful to 
extract a trivial factor of i in order to avoid complex 
matrix elements. Alternatively, by using the relations 



ha 



(at) 
a) a + aa) 



(34) 



valid in the hard-core limit, one can directly decouple the 
kinetic energy term in (|3l|): 



Arta a . 



= 1 
= 1 



Art a] a. 
At 



cosh 



i(«I+«,) 2 
2 ^ 

cr=±l 



(35) 



The same decoupling is used for the hermitian conjugated 
term aja^ yielding, together with the SfS* decoupling, 
three auxiliary Ising-type fields for the XXZ model which 
can be treated in analogy to the algorithm described in 
Sec. II. 



B. Results 



quadratically with At to the exact value, however with 
different prefactors. 
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FIG. 3. AF Heisenberg model on two sites. Full enumer- 
ation over all possible auxiliary field configurations. At ex- 
trapolation with 2nd order Trotter decomposition. 
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□ XY, At=1/4, decoupling (6) 
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FIG. 4. Average sign of the weight function for AF Heisen- 
berg and XY models on four lattice sites (J = 1,B = 0). a 
decreases exponentially with inverse temperature (3. 

Even for very small systems at relatively high temper- 
ature severe minus-sign problems occur. Fig. H shows the 
average sign (a) of w vs. (3 for systems with four sites. 
(a) vanishes exponentially with /3. Values of (a) below 
approx. 0.2 preclude an efficient Monte Carlo sampling. 
The values of (a) do not much differ for the two different 
decouplings, (j?]) and (|35|), of the XY contribution. 

Hence the algorithm does not appear to be suitable for 
models with couplings in more than one spin component. 



In order to check the algorithm and the code we cal- 
culate the nearest-neighbor spin-spin correlation function 
of the isotropic antiferromagnetic (AF) Heisenberg model 
on a 2 x 2 lattice. First, we did a full enumeration over 
all possible auxiliary spin configuration for several val- 
ues of At. Fig. |^ shows the transverse and longitudinal 
correlation function vs. At 2 . Apparently they converge 



IV. CONCLUSIONS 

We have formulated an auxiliary field Quantum Monte 
Carlo algorithm for spin and hard-core boson systems. In 
boson language, it is based on a Hubbard-Stratonovich 
decoupling of the kinetic energy term, leaving a set of 



G 



independent one-site problems in a fluctuating external 
field. Such a procedure has been used in analytic studies 
of the boson-Hubbard model || . The algorithm scales 
linearly with the spatial lattice size, and inverse temper- 
ature, sharing that attractive feature of world-line ap- 
proaches compared to fermion auxiliary field techniques. 
However, unlike traditional world-line techniques it has 
very short auto-correlation times. 

Unfortunately, like fermion auxiliary field approaches, 
the determinants can go negative, resulting in a sign 
problem. We showed in the case of the Ising model in 
a transverse field that an appropriate spin rotation can 
eliminate the problem, making our approach a valuable 
one for studying this problem in more that (1+1) di- 
mensions, where the Jordan- Wigner approach does not 
work. Traditionally formulated world-line simulations, 
which would map the problem onto a highly anisotropic 
classical Ising model, suffer from large auto-correlation 
times that are absent in the present algorithm. Finally, 
our approach shows promise for the extraction of dynami- 
cal correlation functions, a key problem in understanding 
glassy dynamics in the random field case. 
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